#### Meeting-level analysis that propagates the prediction uncertainty ####

source("./code/loadPackages.R")      # Load and install necessary packages
source("./code/regressIter.R")       # Code to run a separate regression for each bootstrap iteration and aggregate results
source("./code/makeCoefPlotsIter.R") # Code to create appropriate coefficient plots for this analysis

# Load in the appropriate meeting-level data
nscdata = fread("./data/meetingLongData.csv")

# Define control variables
controls = " + nAttend + numDef + numInt + numMil + numState + 
              log(totDip+1) + log(totInt+1) + log(totMil+1) +  
              log(usmidschallenge5+1) + cinc + formal +
              mention.Economy + mention.Europe + mention.Asia + 
              mention.Defense + mention.China + mention.Other + 
              mention.Intelligence + mention.Middle.East + mention.Strategic.Forces + 
              mention.Organization + mention.USSR + mention.Policy + 
              mention.Americas + mention.International.Institutions + mention.Diplomacy + 
              mention.Vietnam + mention.Africa + mention.Arms.Control + 
              mention.North.Africa + mention.Latin.America + factor(admin)"
controlsNoAdmin = gsub("\\+ factor\\(admin\\)", "", controls)

# A function to help resolve some errors after running parallel processes
unregister_dopar <- function() {
  env <- foreach:::.foreachGlobals
  rm(list=ls(name=env), pos=env)
}

# Set up cores for parallel processing
no_cores = detectCores() - 1
cl = makeCluster(no_cores)
registerDoParallel(cl)

# Create difference variable for the OLS regressions
nscdata$nDiffConfCoop = nscdata$nConfAdv - nscdata$nCoopAdv

## Emergence model
meanFit1ap = regressIter(outcome="nConfAdv", iv1="meanHawk", controls="+ formal", theFamily="poisson", dataset=nscdata)
meanFit1ap

meanFit1bp = regressIter(outcome="nConfAdv", iv1="meanHawk", controls=controls, theFamily="poisson", dataset=nscdata)
meanFit1bp

meanFit1cp = regressIter(outcome="nDiffConfCoop", iv1="meanHawk", controls="+ formal", theFamily="gaussian", dataset=nscdata)
meanFit1cp

meanFit1dp = regressIter(outcome="nDiffConfCoop", iv1="meanHawk", controls=controls, theFamily="gaussian", dataset=nscdata)
meanFit1dp


## Leader model
leaderFit1ap = regressIter(outcome="nConfAdv", iv1="presHawk", controls="+ formal", theFamily="poisson", dataset=nscdata)
leaderFit1ap

leaderFit1bp = regressIter(outcome="nConfAdv", iv1="presHawk", controls=controlsNoAdmin, theFamily="poisson", dataset=nscdata)
leaderFit1bp

leaderFit1cp = regressIter(outcome="nDiffConfCoop", iv1="presHawk", controls="+ formal", theFamily="gaussian", dataset=nscdata)
leaderFit1cp

leaderFit1dp = regressIter(outcome="nDiffConfCoop", iv1="presHawk", controls=controlsNoAdmin, theFamily="gaussian", dataset=nscdata)
leaderFit1dp


## Adviser model 
compFit1ap = regressIter(outcome="nConfAdv", iv1="advHawk", iv2="presHawk", controls="+ formal", theFamily="poisson", dataset=nscdata)
compFit1ap

compFit1bp = regressIter(outcome="nConfAdv", iv1="advHawk", iv2="presHawk", controls=controlsNoAdmin, theFamily="poisson", dataset=nscdata)
compFit1bp

compFit1cp = regressIter(outcome="nDiffConfCoop", iv1="advHawk", iv2="presHawk", controls="+ formal", theFamily="gaussian", dataset=nscdata)
compFit1cp

compFit1dp = regressIter(outcome="nDiffConfCoop", iv1="advHawk", iv2="presHawk", controls=controlsNoAdmin, theFamily="gaussian", dataset=nscdata)
compFit1dp


## Adviser model with administration fixed effects, without leader hawkishness
compFit2ap = regressIter(outcome="nConfAdv", iv1="advHawk", controls="+ formal + factor(admin)", theFamily="poisson", dataset=nscdata)
compFit2ap

compFit2bp = regressIter(outcome="nConfAdv", iv1="advHawk", controls=controls, theFamily="poisson", dataset=nscdata)
compFit2bp

compFit2cp = regressIter(outcome="nDiffConfCoop", iv1="advHawk", controls="+ formal + factor(admin)", theFamily="gaussian", dataset=nscdata)
compFit2cp

compFit2dp = regressIter(outcome="nDiffConfCoop", iv1="advHawk", controls=controls, theFamily="gaussian", dataset=nscdata)
compFit2dp

# Stop the clusters
stopCluster(cl)
unregister_dopar()


## One single table with only full models
makeFullCoefTable(list(meanFit1bp, meanFit1dp, leaderFit1bp, leaderFit1dp,
                       compFit1bp, compFit1dp, compFit2bp, compFit2dp))

#### Table A15: Effect of Mean Participant Hawkishness and President's Hawkishness on Foreign Policy Decisions, Propagating Uncertainty from Bootstrapping ####
makeFullCoefTable(list(meanFit1ap, meanFit1bp, meanFit1cp, meanFit1dp, 
                       leaderFit1ap, leaderFit1bp, leaderFit1cp, leaderFit1dp))

#### Table A16: Effect of Adviser Hawkishness on Foreign Policy Decisions, Propagating Uncertainty from Bootstrapping ####
makeFullCoefTable(list(compFit1ap, compFit1bp, compFit1cp, compFit1dp,
                       compFit2ap, compFit2bp, compFit2cp, compFit2dp))

#### Figure A9: Summary of Three Models of Trait Aggregation, Propagating Uncertainty from Bootstrapping ####
plotHawkCoefs2("p")
ggsave(filename="./figures/coefPlotUncertain.pdf", height=4, width=7, units='in')
plotMainHawkCoefs2("p")
ggsave(filename="./figures/coefPlotMainUncertain.pdf", height=4, width=7, units='in')
